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ABSTRACT 


The  satellite  kinematics  model  in  the  Naval  Wargaming  System 
(NWGS)  is  described  and  critiqued.  Suggestions  are  made  for 
simplifying  the  existing  computer  code  and  for  increasing  the 
accuracy  of  the  satellite  simulation. 


A.  Introduction 

The  M10  satellite  control  model  (SCM)  in  the  Naval  Warfare 
Gaming  System  (NWGS)  serves  two  primary  purposes.  First,  the  model 
performs  orbital  kinematic  computations  to  position  the 
satellite (s)  to  the  current  game  time.  Second,  the  model  selects 
either  the  radar  satellite  detection  model  or  the  elint  detection 
model  depending  upon  the  type  of  satellite  being  processed.  The 
purpose  of  this  report  is  to  critique  the  orbital  kinematics 
section  of  the  SCM  and  to  make  recommendations  for  improvement. 

B.  The  SCM  Kinematics 

Satellite  parameter  inputs  to  the  SCM,  pertinent  to  the 
kinematics,  consist  of  the  satellite  orbital  period,  the  time  of 
last  position  update  and  the  latitude,  longitude  and  course  at  the 
time  of  last  position  update.  The  time  interval  since  the  last 
position  update,  t_int,  is  the  current  game  time  minus  the  time  of 
last  update.  The  great  circle  distance  traveled  (measured  in  arc 
length)  since  the  time  of  last  update  is  t_int/period.  Using 
spherical  trigonometry,  the  position  and  course  are  updated  by 
moving  the  satellite  along  the  great  circle  determined  by  its  last 
position  and  course,  and  the  arc  distance  traveled  in  t_int.  The 
ground  position  is  then  corrected  by  subtracting  the  amount  of  the 
Earth's  rotation  in  t_int  from  the  updated  longitude. 

In  essence,  the  SCM  assumes  all  satellites  to  be  in 
unperturbed  circular  Keplerian  orbits  about  the  Earth. 
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C.  Critique 

The  major  limitation  of  the  SCM  is  that  all  orbits  are  assumed 
circular.  Even  when  this  assumption  is  valid,  the  major  flaw  in 
the  existing  SCM  is  the  neglect  of  the  precession  of  the  orbit 
caused  by  the  Earth's  oblateness.  This  precession  depends  upon  the 
height  and  inclination  of  the  orbit  and  can  be  as  much  as  10 
degrees  per  day  which,  in  turn,  can  cause  a  ground  position  error 
up  to  600  n.  mi.  per  day.  Neither  orbital  height  nor  inclination 
are  input  parameters  to  the  SCM. 

An  additional  error,  caused  by  the  neglect  of  the  Earth's 
oblateness,  is  in  the  computed  latitude  of  the  satellite's  ground 
track.  At  current  game  time  the  computed  (geocentric)  latitude  can 
be  in  error  by  as  much  as  11.5  n.  mi.  from  the  true  (geographic) 
latitude  at  latitudes  near  45  degrees. 

D.  Summary  of  Recommendations 

The  sophistication  of  the  SCM  model  can  be  increased  in 
several  steps.  An  overview  is  given  below.  Details  are  given  in 
Section  E. 

1.  The  programming  of  the  existing  SCM  can  be  simplified  by 
using  alternate  formulas  from  spherical  trigonometry. 

2.  By  including  the  satellite  semimajor  axis  (which  can  be 
derived  from  its  orbital  period)  and  the  orbital 
inclination,  the  existing  SCM  can  be  easily  modified  to 
include  the  orbital  precession  caused  by  the  Earth's 
oblateness.  Without  this  inclusion,  a  ground  position  error 
of  600  n.mi.  per  day  can  occur. 
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3.  To  include  noncircular  orbits  it  will  be  necessary  to 
introduce  all  of  the  satellite  orbital  parameters.  The 
computations  involve  the  solution  of  Kepler's  equation,  the 
transformation  of  the  in-plane  coordinates  to  an  inertial 
system,  and  the  reduction  of  the  coordinates  to  geocentric 
latitude  and  longitude. 

4.  The  satellite's  geocentric  latitude  can  be  converted  to 
geographic  latitude  to  account  for  the  Earth's  oblateness. 
This  conversion  can  eliminate  an  11.5  n.mi.  ground  position 
error  in  the  N-S  direction.  A  by-product  is  the 
determination  of  the  satellite's  height. 

The  computation  of  the  satellite  position  using  Recommenda- 
tions 3  and  4  are  time  consuming  and  could  degrade  real  time 
operation  of  the  existing  NWGS.  It  would,  however,  be  possible  to 
precompute  and  store  the  satellite  positions  in  a  table  prior  to 
running  the  NWGS. 

E.  Detailed  Recommendations 

1.     Add  a  two  argument  quadrant  arctangent  function,  called 
here  'qatanbf,  to  the  atrigf  include  file.  Details  of  the 
qatan  function  are  in  the  appendix.  Then,  the  Ml0_satellite_ 
control_ee  module  can  be  modified  as  follows: 
a.  Replace  lines  201-202  by 

d_long  =  qatanbf  (sinb_dist  *  sinb„course  ,  cosb_dist  * 
cosb_lat  -  sinb_dist  *  cosb_course  *  sinb„lat) . 
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b.  Replace  lines  206-220  by 

n_course  =  qatanbf  (  sinb_course  *  cosb_lat  , 

cosb_course  *  cosb_dist  -  sinb_dist  *  sinb_lat) . 

c.  With  these  two  modifications,  lines  224-231  are  no  longer 
required. 

These  equations  take  advantage  of  the  general  spherical 
triangle  [Ref .  1] .  Both  d_long  and  n_course,  as  computed  above, 
lie  in  the  range  of  -1  to  +1  bams  (-pi  to  +pi) .  The  4-part 
formula  [Ref.  2]  is  used  to  compute  n_course  directly.  The 
possibility  of  polar  crossing  is  automatically  included. 

2.     The  formula  for  the  precession  of  the  ascending  node  of  the 
orbit,  owing  to  the  Earth's  oblateness,  is 


(3  J9  cos  i)n(T  -  Ta) 

N  N0  ~  2 

4.Q. 


where  N  is  the  longitude  of  the  ascending  node  at  time  T,  N,,  is 
the  longitude  of  the  ascending  node  at  time  T.,,  J~  = 
0.001082616  is  the  Earth's  second  harmonic  shape  parameter,   i 
is  the  orbital  inclination,   n   is  the  satellite  motion  and  a 
is  the  semimajor  axis.  The  motion   n  =  2  /P,  where  P  is  the 
orbital  period.  Further,   a   and   P   are  related  by 

«2      .  2  3.2 
P   =   4tt  a  /k  , 

where  k  =  0. 0743669161 (e. r. ) 3'2/min.  is  the  Earth's 
gravitational  constant.  Using  these  relations,  the  formula  for 
the  node  may  be  written  in  one  of  several  ways: 
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N  =  N0  -  9.9639(ae/a)7/2(T  -  T0)  cos  i, 

where  a   is  the  Earth's  equatorial  radius,  N  is  degrees  and  T 
is  measured  in  days,  or 

N  =  N0  -  216.74  P~7/3(T  -  T0)  cos  i, 

where  N  is  in  degrees,  and  P  and  T  are  measured  in  minutes,  or, 

N  =  N0  -  282.83  P~7//3t_int  cos  i, 

where  N  is  in  bams,  and  P  and  t_int  are  in  seconds  as  is  the 
custom  in  the  NWGS. 

Using  the  last  equation,  line  235  of  Ml0_satellite_ 
control_ee  can  be  written  as 

sat_long  =  sat_long  +  d_long  -  (bams_sec  +  prec)  *  t_int, 

-7/3 
where   prec  =  282.83  *  P  '"  *   cos  i.  Note  that   prec   is  a 

constant  for  each  satellite. 

3.     Since  noncircular  orbits  do  not  give  rise  to  uniform  motion 
of  the  ground  track  over  a  non-rotating  spherical  Earth,  the 
method  used  in  M10_satellite_control_ee  cannot  be  used.  In 
fact,  eccentric  orbits  can  give  rise  to  ground  tracks 
containing  cusps. 

Whatever  procedure  is  used  for  the  orbital  computation,  it 
will  be  necessary  to  obtain  all  of  the  satellite  parameters. 
These  are:  the  semimajor  axis   a,  the  eccentricity   e,  the 
longitude  of  the  ascending  node   N0,  the  argument  of  perigee 
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w0,  the  angle  of  inclination   i,  and  the  mean  anomaly  M0   for 
some  epoch  T0,  or  their  equivalent.  Because  of  the  complexity 
of  the  computations,  it  might  be  advisable  to  obtain  the 
orbital  elements  distributed  by  NORAD  and  to  use  simplest  of 
their  existing  programs,  the  SGP  model  [Ref .  3] .  The  NORAD 
elements  differ  slightly  from  those  given  above  but  they  also 
compensate  for  orbital  decay  due  to  atmospheric  drag.  Another 
alternative  would  be  to  use  equations  such  as  those  given  in 
Escobal  [Ref.  4].  An  outline  of  the  necessary  computations  is 
given  below. 

Set  up  the  global  constants: 

k  =  0.0743669161  =  gravitational  constant, 
J2  =  0.001082616  =  second  harmonic  of  Earth, 
a  =  1  e.r.  =  Earth's  equatorial  radius, 
T'  =  0.0043752695  rad/min  =  Earth's  rotation  rate. 

Input  a,   e,   i,   w0,   N0,   M0,  and  T0. 
Compute  the  satellite  constants: 

p  =  a(l  -  e2), 

n  =  ka"3/2[l  +  1.5  J2(l  -  e2)"1/2(l  -  1.5  sin2  i)/p2] , 

w*  =  1.5  J2  n(2  -  2.5  sin2  i)/p2, 

N'  =  -  (1.5  J2  n  cos  i)/p2. 

The  time  interval,  measured  in  minutes  from  T0,  is  T  -  T0, 
where  T  is  the  current  time.  Compute: 
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w  =  w0  +  w" (T  -  T0) , 


N  =  N  +  N' (T  -  T0) , 


M  =  mod[n(T  -  T0)  ,  2 tt ]  . 


Solve  Kepler's  equation  for  E 

E  -  e  sin  E  =  M. 


Compute  the  in-plane  position  components,  X   and  Y  ,  and 

W  W 

the  radius  vector  R: 


X  =  a (cos  E  -  e) , 
w 

Y  =  a(l  -  e2)1/2  sin  E, 
w 

R  =  a (1  -  e  cos  E) . 
Reduce  the  in-plane  coordinates  to  equatorial  coordinates 


P   =  cos  w  cos  N  -  sin  w  sin  N  cos  i, 
x 

P  =  cos  w  sin  N  +  sin  w  cos  N  cos  i, 

y 

P   =  sin  w  sin  i, 

Q   =  -  sin  w  cos  N  -  cos  w  sin  N  cos  i, 
x 

Q  =  -  sin  w  sin  N  +  cos  w  cos  N  cos  i, 

Q   =  cos  w  sin  i, 
z 

X  =  P   X   +  Q   Y  , 
X   w     x   w' 

Y  =  P   X   +  Q   Y  , 

Z  =  P   X   +  Q   Y  . 
z   w     z   w 
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Compute  the  right  ascension,  RA,  and  declination,  DEC,  of 
the  satellite: 

RA  =  qatan(Y,X) , 
DEC  =  arccos(Z/R) . 

The  geocentric  latitude  of  the  satellite  L'  =  DEC. 
Compute  the  unnormalized  longitude: 

L1  =  RA  -  T   -  T't, 
n         g    g 

where   t   is  the  total  number  of  minutes  elapsed  since  the 
year,  month  and  day  of  T0,  and  T   is  the  Greenwich  sidereal 
time  at  the  initial  time  T0.   T  ,  measured  in  degrees,  is  given 
by 

T   =  99.6909833  +  (36000.7689  +  0.00038708T  )T  , 
g  u'  u' 

where  T   is  given  by 

Tu  =  (JD  -  2415020. 0)/36525, 

and  JD  is  the  Julian  Day  Number, 

JD  =  367Y  -  int{7[Y+int( (N  +  9)/12)]/4} 

+  int(275M/9)  +  D  +  1721013.5  +  UT/24, 

where  Y  is  the  year  (1901  1  Y  i.  2099),  M  is  the  month  (1  X  M  £ 
12),  D  is  the  day  (1  1  D  <,   31),  and  UT  is  the  universal  time  in 
hours.  Then  compute  the  east  longitude  L,  where 


Ln  =  mod(L' ,  2tt)  . 
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4.     A  conversion  of  geocentric  latitude,  L",  to  geographic 

latitude  for  an  oblate  Earth  can  be  found  in  Ref .  4.  A  slightly 
improved  version  is  given  below.  Let 

a   =  6378.135  n.  mi.  =  Earth's  equatorial  radius, 

f  =  1/298.26  =  Earth's  flattening  factor, 

r  =  a  R. 
e 

Let  L'  =  DEL.  Then  compute 

rc  =  ae{[l  -  (2f  -  f2)]/U  -  (2f  -  f2)  cos2  L']}1/2' 
Lfc  =  arctan[tan  L'/(l  -  f)2], 


H/r  =  [1  -  (rc/r)2  sin2(Lt  -  L']1/2 


-  (rc/r)  cos(Lt  -  L') , 

DL'  =  arcsin[(H/r)  sin(Lfc  -  L')j. 

Next,  let  L'  =  DEL  -  DL '  and  recompute  from  the  equation  for  r 
until  L'  no  longer  varies  (convergence  is  rapid).  L.  is  the 
geographic  longitude  of  the  satellite  and  H  =  r(H/r)  is  the 
height  of  the  satellite  above  the  oblate  Earth. 


References 

1.  Chauvenet,  William,  "Treatise  on  Plane  and  Spherical 
Trigonometry" ,  1850. 

2.  Smart,  W.  M.,  "Spherical  Astronomy",  Cambridge  University 
Press,  1965. 

3.  Hoots,  Felix  R.  and  Roehrich,  Ronald  L.,  "Models  for 
Propagation  of  NORAD  Element  Sets",  Spacetrack  Report  No.  3, 
Project  Space  Track,  Aerospace  Defense  Command,  United  States 
Air  Force,  December  1980. 

4.  Escobal,  Pedro  R. ,  "Methods  of  Orbit  Determination",  John  Wiley 
&  Sons,  1965. 


-10- 


Appendix 

The  quadrant  arctangent  function,  designated  qatan,  is  a 
two  argument  function  of  rectangular  coordinates   x   and  y 
which  returns  an  arctangent  value  in  the  range  of  -it  to  +it 
depending  upon  the  quadrant  of   x   and   y.  In  FORTRAN  this 
function  is  called  ATAN2.  It  can  be  computed  as  follows: 


p  ~.  c 


qatan  (x,y)  =  atan{y/[x  +  e*t  (x— ^-0-)  J  } 

+  7T*t(x  <  0)*[sign(y)  +  t(y  £   0)  ]  , 


where 


e  =  smallest  number  representable  on  the  computer, 

t(z)  =  1  when  z  is  true, 
t(z)  =  0  when  z  is  false, 


sign(y)  = 


+1  if  y  >  3, 

3  if  y  =  0, 

-1  if  y  <  0. 
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